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A new numerical method for the solution of the Dynamical Mean Field Theory's self-consistent 
equations is introduced. The method uses the Density Matrix Renormalization Group technique 
to solve the associated impurity problem. The new algorithm makes no a priori approximations 
and is only limited by the number of sites that can be considered. We obtain accurate estimates 
of the critical values of the metal-insulator transitions and provide evidence of substructure in the 
Hubbard bands of the correlated metal. With this algorithm, more complex models having a larger 
number of degrees of freedom can be considered and finite-size effects can be minimized. 
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Great theoretical progress in our understanding of the 
physics of strongly correlated electron systems has been 
possible since the introduction of the Dynamical Mean 
Field Theory (DMFT) just over ten years now [ll|3|. This 
approach is based on the natural extension of the famil- 
iar classical mean-field theory of statistical mechanics to 
the treatment of models of strongly interacting electrons 
on a lattice. The DMFT solution of the model is exact 
in the limit of large lattice dimensionality or large con- 
nectivity 2, 3]. Since its introduction, DMFT has been 
widely adopted and was used for the investigation of a 
large variety of model Hamiltonians relevant for problems 
as diverse as colossal magneto-resistance, heavy fermions, 
metal-insulator transitions, etc. A great deal of in- 
terest is currently centered around the ongoing efforts 
to incorporate DMFT as the local correlation physics 
"engine" for first-principle calculations of realistic com- 
pounds nil]. At the heart of the DMFT method is the 
solution of an associated quantum impurity model where 
the environment of the impurity has to be determined 
self-consistently. Therefore the ability to obtain reliable 
DMFT solutions of lattice model Hamiltonians relies di- 
rectly on the ability to solve quantum impurity models. 
Since solutions of general impurity models are usually 
not analytically tractable, one has to resort to numeri- 
cal algorithms or approximate methods. Among the a 
priori exact numerical algorithms available we count the 
Hirsch-Fye Quantum Monte Carlo 6] method and Wil- 
son's Numerical Renormalization Group (NRG) • The 
former is a finite-temperature method that is formulated 
in imaginary time and has been applied to a large variety 
of impurity models including the multi-orbital case that 
corresponds to correlated multi-band lattice Hamiltoni- 
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ans While this method is very stable and accurate 
having allowed for extremely detailed investigations of 
fundamental problems such as the metal-insulator tran- 
sition in the Hubbard model , its main drawback is 
that the access to real frequency quantities such as spec- 
tral functions requires to recourse to less controlled tech- 
niques for the analytic continuation of the Green func- 
tions. The second numerical method is based on Wilson's 
renormalization group [lOl . This method can be for- 
mulated both at T=0 and finite (small) T and provides 
extremely accurate results at very small frequencies, thus 
it is very adequate for the investigation of correlated 
metallic phases with heavy effective-mass quasi-particles. 
The cost to pay is, however, that the description of the 
high energy features involve some degree of approxima- 
tion and cannot be so accurately obtained [l^ . 

The goal of the present work is to introduce a new 
algorithm for the solution of the DMFT self-consistent 
equations, that makes use of another powerful numeri- 
cal methodology for the solution of many-body Hamil- 
tonians: the Density Matrix Renormalization Group 
(DMRG)[i3|. This method, fike the NRG, has the ap- 
pealing feature of making no a priori approximations and 
the possibility of a systematic improvement of the quality 
of the solutions. However, unlike NRG, it is not formu- 
lated as a low-frequency asymptotic method and thus 
provides equally reliable solutions for both gapless and 
gapfuU phases. More significantly, it provides accurate 
estimates for the distributions of spectral intensities of 
high frequency features such as the Hubbard bands, that 
are of main relevance for analysis of x-ray photoemission 
and optical conductivity experiments. 

We shall illustrate the new formulation with the solu- 
tion of the, by now classic, Mott transition in the Hub- 
bard model. We shall show that accurate estimates of the 
critical values of the interaction for the metal-to-insulator 
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and for the insulator-to-metal transitions at T=0 can be 
obtained, and, more interestingly, we shall provide evi- 
dence of substructure in the Hubbard bands in the cor- 
related metallic phase. 

The Hamiltonian of the Hubbard model is defined by 

H = tJ2 4<tCj\<t + f^Xl'^^i"*.T (1) 

The treatment of this model Hamiltonian with DMFT 
leads to a mapping of the original lattice model onto an 
associated quantum impurity problem in a self-consistent 
bath. In the particular case of the Hubbard model, 
the associated impurity problem is the single impurity 
Anderson model (SIAM), where the hybridization func- 
tion A (a;), which in the usual SIAM is a flat density of 
states of the conduction electrons, is now to be deter- 
mined self-consistently. More precisely, for the Hamilto- 
nian defined on a Bethe lattice of coordination d, one 
takes the limit of large d and exactly maps the model 
onto a SIAM impurity problem with the requirement 
that A{ui) = t^G{ijj), where G{uo) is the impurity Green's 
function. At the self-consistent point G{lo) coincides with 
the local Green's function of the original lattice model 0. 
We take the half-bandwidth of the non interacting model 
as unit of energy, thus t = 1/2. 

A central quantity in this algorithm is the non in- 
teracting Green's function of the impurity problem, 
Go(u;) = l/(w + iJL- A(w)) = l/(u; + - t^G[Lo)). 
Thus, to implement the new algorithm we shall con- 
sider 0, 01 a general representation of the hybridiza- 
tion function in terms of continued fractions that de- 
fine a parametrization of A{lo) in terms of a set of 
real and positive coefhcients. Since it is essentially a 
Green's function, A(z) can be decomposed into "parti- 
cle" and "hole" contributions as A{z) = A>(z) -I- A<(z) 
with A>(z) = t^{gs\cj^^j^^c^\gs) and A<{z) = 

t'^{gs\c^ c\gs) for a given Hamiltonian, H with 

ground-state energy Eq. By standard Lanczos technique, 
H can be in principle tri-diagonalized and the functions 
A>(z) and A<(z) can be expressed in terms of respec- 
tive continued fractions^]. As first implemented in 
Ref.[ll[li, each continued fraction can be represented 
by a chain of auxiliary atomic sites whose energies and 
hopping amplitudes are given by the continued fraction 
diagonal and off-diagonal coefficients respectively. 

From the self-consistency condition, the two chains 
representing the hybridization, are "attached" to the 
right and left of an atomic site to obtain a new SIAM 
Hamiltonian, H. In fact Go{z) constitutes the local 
Green's function of the site plus chain system. 

The algorithm in Ref.0, basically consists in 
switching on the local Coulomb interaction at the impu- 
rity site of the SIAM Hamiltonian and use the Lanczos 
technique to re-obtain A(2;), iterating the procedure until 
the set of continued fractions coefhcients converges. 



A great limitation of this procedure is that the num- 
ber of auxiliary atomic sites that needs to be used in 
the hybridization chains is too large for standard exact 
diagonalization schemes. Therefore, the chains have to 
be cut at very short lengths and as a consequence the 
method becomes approximate with sizable finite-size ef- 
fects. In the present work, we shall make a fundamental 
improvement on this scheme, namely diagonalizing the 
Hamiltonian using DMRG, which in principle allows to 
handle chains of arbitrary length [l3| . 

The SIAM Hamiltonian therefore reads 

Nc Nc-1 

H ^ ^ aaC^^„Caa+ ^ ba{cl^„Ca+la + h.C.) 

<7,a=~Nc cr.a = — {Nc — l) 

a=£0,-l 

+ bo{cicc.„ + h.c.) + U{n^~^){ni^^) (2) 

(7,a— ±1 

with Ca being the destruction operator at the impurity 
site, and Caa being the destruction operator at the a site 
of the hybridization chain of 2Nc sites. The set of param- 
eters {aa,ba} are directly obtained from the coefficients 
of the continued fraction representations of A(z) by the 
procedure just described. 

An important point to make is that while the hopping 
to an impurity connected to a chain bears a resemblance 
to the NRG method, the present formulation is different 
in a crucial aspect: unlike NRG, it is not constructed as 
an asymptotically exact representation for low frequen- 
cies, but rather treats all energy scales on equal footing. 
By paying the price of giving up the excellent resolution 
at low frequency of NRG, we shall gain in exchange a 
controlled and systematic algorithm operating at all en- 
ergy scales. Of course, as in Wilson's NRG the energy 
resolution depends on the length of the auxiliary chain, 
therefore considering longer chains while keeping the nu- 
merical accuracy in the computation is the central limi- 
tation of the current scheme. In practice, systems of up 
to 45 {Nc = 22) sites have been considered. 

We now turn to the results of the implementation 
of this algorithm for the solution of the DMFT equa- 
tions of the Hubbard model. In FigQ] we show the 
DMFT+DMRG results (solid fines) for the density of 
states (DOS) for several values of increasing interaction 
U. The results are compared to the Iterated Perturba- 
tion Theory (IPT) results (dashed fines) 0, 13 which is 
a useful analytic approximate method that can be solved 
on the real frequency axis at T = 0. In Fig|21we compare 
the DMFT-I-DMRG results (solid lines) for the imaginary 
part of the Green's function on the imaginary frequency 
axis with the corresponding ones obtained by quantum 
Monte Carlo solution of the DMFT equations at a low 
temperature (circles). We see that the overall agreement 
on the real axis is very satisfactory and on the Matsubara 
frequency axis it is excellent. At low values of the interac- 
tion U we find a metallic state characterized by a narrow 
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quasi-particle feature at low frequencies. Increasing U 
this peak gets narrower by transferring spectral weight 
to features at higher energies of order U, the upper and 
lower Hubbard bands. At large values of the interaction 
the system evolves toward an insulating state with a gap 
of order U (Fig. I^c)). 
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FIG. 1: Density of states {^lmG{uj)) for the half-filled Hub- 
bard Model ,18^] . We also show the IPT results (dashed lines). 
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FIG. 2: Imaginary part of the Green functions at imaginary 
(Matsubara) frequencies (solid lines) . We also show Quantum 
Monte Carlo results (circles) at low temperature T — 1/32. 

A very important feature of the metal-insulator transi- 
tion in the paramagnetic state of the Hubbard model at 
half-filling (2i is that there are two distinct critical values 
of the interaction associated with the transition: Ud and 
Uc2- The former signals the insulator to metal transition 
obtained upon lowering the interaction, while the latter 
corresponds to the metal to insulator transition obtained 
when the Fermi liquid is destroyed by increasing the in- 
teraction strength. We obtained estimates of these two 
values that are consistent with those from NRG calcula- 
tions [13. We find Uci = 2.39 ±0.02 and Uc2 = 3.0 ±0.2. 
Due to the nature of this algorithm and the arguments 



presented before, we expect that our determination of 
Uci should be more accurate than NRG (and all other 
previously used T = methods). Our criterion for the in- 
vestigation of metallic versus insulating states was based 
on the behavior of the spectral weight at zero frequency 
and the size of the gap in the DOS (given by the en- 
ergy of the first pole). It was a remarkable finding that 
these quantities showed an unexpected dramatic depen- 
dence with the length of the hybridization chain and the 
proximity to the critical value of the interaction. This 
dramatic effect is demonstrated in Fig|3| where we plot 
the results as a function of the inverse of the chain length 
at [/ = 2.3 (a) and U = 2.5 (b). We find that as we in- 
crease the value of U from the weak coupling (metallic) 
side, chains longer than a [/-dependent value were re- 
quired to converge to metallic solutions (Fig|3|d)). The 
value of Lc was in fact found to diverge at a finite in- 
teraction strength which we identified with our estimate 
of Uc2- As long as the length L = Nc < Lc the solu- 
tion looks as that of an insulating state with vanishing 
density of states at lj = 0, while a rapid crossover to a 
metallic solution is seen as L goes beyond the threshold 
value Lc- Interestingly, this resembles the behaviour of 
the Kondo effect in finite systems ,20|, where the Kondo 
effect is strongly suppressed if the Kondo screening cloud 
is larger than the system size. This similarity also shows 
up in the increase of Lc with U , since in the Kondo model 
the correlation length increases with the interaction. 

On the other hand, for the insulator to metal transi- 
tion, the value of Uci can be determined by the closure 
of the gap in the DOS or using the inverse of the second 
moment of the DOsHJ], Fig. Elc). Any of these two 
quantities are non zero in the insulating state and vanish 
at Uci- 

Another interesting result that has been a matter of 
debate, and might have implications for the analysis of 
x-ray photoemission spectra, is the question of the exis- 
tence on substructure in the Hubbard bands of the cor- 
related metallic state. The substructure was first identi- 
fied within IPT calculations but the exact numerical 
methods did not have the required accuracy to decide 
whether the substructure was an artifact borne out of 
IPT or a real feature of the model's solution. In fact, an 
appealing physical interpretation of the substructure can 
be readily made: in a rigid band picture, one can think of 
the action of the local interaction U to "split" and repli- 
cate the hybridization function density of states A{lu) at 
frequencies ±C//2. This in fact gives a simple qualitative 
understanding of the emergence of Hubbard bands in the 
insulator when U is large. The semicircular A(a;) is du- 
plicated into the two (approximate) semicircles at ±[//2 
of the local DOS (the imaginary part of G(uj)) character- 
istic of the Mott-Hubbard insulator. Then one realizes 
that the self-consistency requires A{uj) to coincide with 
t^G{uo), thus, if Im[G(w)] develops (multi-peak) structure 
it implies that A(ti;) has a similar structure. Then, by 
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FIG. 3: Density of states at the Fermi energy, ilmG(cj = 0), 
(squares) and gap (circles) for a) U = 2.3 and h) U = 2.5. 
Open symbols correspond to metallic solutions and closed 
symbols in (b) to the insulating phase. The thick arrow shows 
the critical length (Lc) above which a metallic solution is ob- 
tained, c) Gap (circles) and inverse of the second moment 
of the DOS (diamond) as a function of interaction U . The 
values correspond to extrapolations to infinite-size chains, d) 
Critical length Lc as a function of interaction U . 



the simple-minded argument in which the action of U is 
to "split and replicate" , we get (multi-peak) structures 
for the Hubbard bands at ztU/2 in the DOS, in a self- 
consistent manner. In Fig^we show the comparison of 
the DOS of the model for both the correlated metallic 
and insulating solutions at the same value of the inter- 
action U=2.5 (i.e., between Ud and Uc2)- The results 
clearly show that the smooth envelope of the pole struc- 
ture (due to a finite Nc) forming the lower and upper 
Hubbard bands in the insulator, acquires more structure 
in the metallic state. For further comparison, we also 
plot the corresponding results from NRG|2^ and IPX. 
The former fails to capture any qualitative difference 
between the shape of the metallic and insulating Hub- 
bard bands, however, IPT shows a very smooth shape 
for the insulator (cf Fig^-d) and a multi-peak structure 
in the metallic Hubbard bands. The observed structure 
of IPT is in qualitative agreement with the spectral dis- 
tribution of weight born-out from the DMRG calculation. 
To our knowledge this is the first strong evidence of the 
existence of non-trivial structure in the Hubbard bands 
within DMFT. 

To conclude, we have presented a new algorithm to 
solve the DMFT equations of strongly correlated mod- 
els exploiting the DMRG methodology. Large systems 
can be considered and accurate values of the critical in- 
teractions are obtained in agreement with NRG predic- 
tions allowing for a non-trivial test of the accuracy of 
this method. In contrast with NRG, however, this new 
algorithm deals with all energy scales on equal footing 
which allowed us to find interesting substructure in the 
Hubbard bands of the correlated metallic state. This ob- 




FIG. 4: Density of states for U=2.5 (solid line). For compari- 
son NRG (dot-dashed line) and IPT (dashed line) results are 
also shown, a) Insulating solution, b) Metallic solution. 



servation might be useful for the interpretation of high 
resolution photoemission spectroscopies. In addition this 
method can handle more general models having a larger 
number of degrees of freedom and can be generalized to 
finite temperatures. 
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